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Thermal  energy  storage  is  needed  to  improve  the  efficiency  of  solar  thermal  energy  applications  (STEA) 
and  to  eliminate  the  mismatch  between  energy  supply  and  energy  demand.  Among  the  thermal  energy 
storages,  the  latent  heat  thermal  energy  storage  (LHTES)  has  gained  much  attention  because  of  its  high- 
energy  densities  per  unit  mass/volume  at  nearly  constant  temperatures.  This  review  presents  previous 
studies  on  the  numerical  modeling  of  phase  change  materials  (PCMs)  through  a  commercial  computa¬ 
tional  fluid  dynamic  (CFD)  software  and  self-developed  programming  to  study  the  heat  transfer 
phenomena  in  PCMs.  The  CFD  (Fluent)  software  is  successively  used  to  simulate  the  application  of 
PCMs  in  different  engineering  applications,  including  electronic  cooling  technology,  building  thermal 
storage,  and  heating,  ventilation,  air  conditioning  (HVAC).  Using  CFD  software  to  design  LHTES  is 
believed  to  be  an  effective  way  to  save  money  and  time  and  to  deliver  optimization  tools  for  maximum 
efficiency  of  STEAs. 
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1.  Introduction 

Researchers  intensively  studied  the  thermal  energy  storage 
of  PCMs  for  the  last  three  decades  because  of  the  latter’s  high 
thermal  energy  densities  per  unit  volume/mass  and  their  applic¬ 
ability  to  different  engineering  fields  using  wide  temperature 
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ranges.  PCM  thermal  storage  plays  a  key  role  in  improving  energy 
efficiency.  It  limits  the  discrepancy  between  the  energy  supply 
and  the  energy  demand  of  solar  thermal  energy  applications 
(STEAs),  particularly  when  the  STEA  operation  strategy  depends 
solely  on  solar  energy  as  a  main  source.  PCM  thermal  storage 
indicates  high  performance  and  dependability  with  the  advan¬ 
tages  of  high  storage  capacity  and  nearly  constant  thermal  energy 
[i].  Most  STEAs  need  constant  or  near-constant  temperature  for 
high-efficiency  strategies.  Using  latent  heat  thermal  energy  sto¬ 
rage  (LHTES)  as  thermal  energy  storage  can  provide  the  required 
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Nomenclature 

P 

fluid  density  (kg/m3) 

y 

liquid  fraction 

C 

Mushy  zone  constant  (kg/m3  s) 

p 

volumetric  expansion  coefficient  (1/K) 

Cp 

Specific  heat  of  PCM,  (J/kg) 

p 

Dynamic  viscosity  (kg/m  s) 

g 

gravity  acceleration,  (m/s2) 

a 

Phase  volume  fraction 

h 

sensible  enthalpy(J/kg) 

H 

enthalpy  (J/kg) 

Subscripts 

L 

latent  heat  fusion  (J/kg) 

k 

thermal  conductivity  (W/m  K) 

ref 

Reference 

T 

Temperature  (°C  or  K) 

s 

Solidus  of  the  phase  change  materials 

u 

velocity  (m/s) 

1 

Liquidus  of  the  phase  change  material 

Greek  letters 

constant  temperature  that  matches  the  melting  temperature  of 
the  PCMs.  PCMs  are  used  in  different  engineering  fields,  such  as  in 
the  following:  the  thermal  storage  of  building  structures  [3-7]; 
building  equipment,  including  domestic  hot  water;  heating  and 
cooling  systems  [8,9];  electronic  products  [10-12];  drying  tech¬ 
nology  [13],  waste  heat  recovery  [14];  refrigeration  and  cold 
storage  [15-18];  solar  air  collectors[19];  and  solar  cookers  [20]. 
Using  a  CFD  software  to  design  a  LHTES  is  expected  to  be  an 
effective  way  to  save  money  and  time  and  to  deliver  optimization 
tools  for  maximum  efficiency  of  STEAs. 


2.  Numerical  solution  of  PCMs 


The  mathematical  formulation  of  a  phase  transient  known 
as  a  phase  change  or  moving  boundary  is  governed  by  a  partial 
deferential  equation  that  can  be  solved  either  analytically  or 
numerically.  The  analytical  solution  of  PCMs  is  problematic 
because  of  the  nonlinear  phase  front  interfaces,  complex  geome¬ 
tries,  and  nonstandard  boundary  condition;  the  few  analytical 
studies  available  are  on  ID  cases  with  regular  geometries  and 
a  standard  boundary  condition.  Ozisik  [21]  reported  that  the 
numerical  method  for  solving  PCMs  can  be  categorized  as  fixed- 
grid,  variable  grid,  front-fixing,  adaptive  grid  generation,  and 
enthalpy  methods. 

Predicting  the  behavior  of  phase  change  systems  is  difficult 
because  of  its  inherent  non-linear  nature  at  moving  interfaces,  for 
which  the  displacement  rate  is  controlled  by  latent  heat  lost  or 
absorbed  at  the  boundary  [22].  The  heat  transfer  phenomena  in 
solid-liquid  PCMs  can  be  analyzed  using  two  main  methods:  the 
temperature-based  and  enthalpy-based  methods.  In  the  first 
method,  temperature  is  considered  a  sole  dependent  variable. 
The  energy  conservation  equations  for  the  solid  and  liquid  are 
written  separately;  thus,  the  solid-liquid  interface  position  can  be 
tracked  explicitly  to  achieve  an  accurate  solution  for  the  pro¬ 
blems. 


dTs 

dn 


dTi 

dn 


ki  +  pLkVn , 


where  Ts  denotes  the  temperature  in  the  solid  phase,  T\  denotes 
the  temperature  in  the  liquid  phase,  ks  is  the  thermal  conductivity 
of  the  solid  phase,  k\  is  the  thermal  conductivity  of  the  liquid 
phase,  n  is  the  unit  normal  vector  to  the  interface,  and  vn  is  the 
normal  component  of  the  velocity  of  the  interface.  L  is  the  latent 
heat  of  freezing,  as  shown  in  Fig.  1. 

In  the  second  method,  the  solid-liquid  interface  position  need 
not  be  tracked.  Researchers  often  use  the  enthalpy  formulation 
because  of  the  following  advantages:  (1)  the  governing  equations 
are  similar  to  the  single-phase  Eq.;  (2)  no  explicit  conditions  need 


to  be  satisfied  at  the  solid-liquid  interface;  (3)  the  enthalpy 
formulation  involves  the  solution  within  a  mushy  zone,  involving 
both  solid  and  liquid  materials,  between  the  two  standard  phases; 
and  (4)  the  phase  change  problem  can  be  solved  more  easily  [23]. 

^p+V-(pvH)=V-(kAT)+S  (2) 

Where  T  denotes  the  temperature,  k  is  the  thermal  conductivity, 
p  is  the  density  of  the  PCM,  v  is  the  fluid  velocity,  and  H  is  the 
enthalpy,  S  is  the  source  term. 

Dutil  et  al.  [22]  presented  an  intensive  mathematical  and 
numerical  review  of  the  PCM  application  based  on  the  first  and 
second  laws  of  thermodynamics.  Using  252  references,  they 
determined  the  mathematical  and  numerical  methods  applied 
to  solve  heat  transfer  problems  involving  PCMs  for  thermal 
energy  storage,  the  mathematical  fundamentals  of  PCMs,  and 
the  different  application  geometries  and  applications.  Verma  et  al. 
[24]  also  introduced  other  PCM  mathematical  reviews. 

Recently,  researchers  used  the  Fluent  software  by  ANSYS  to 
simulate  melting  and  solidification  in  engineering  problems. 
Other  software  that  can  be  used  to  simulate  the  PCM  process 
include  COMSOL  Multiphysics  and  Star-CMM  +  .  However,  Fluent 
is  preferred  by  most  researchers  for  melting  and  solidifying  PCMs. 
Conversely,  some  researchers  self-developed  a  program  using 
computational  language  (C-H-,  Fortran,  Matlab)  to  study  the  heat 
transfer  phenomena  in  PCMs.  Table  1  summarizes  some  of  the 
self-developed  programs  for  different  PCM  geometries. 

2.2.  Fluent  program 

The  Fluent  software  by  ANSYS  is  a  computational  fluid 
dynamic  (CFD)  program  used  successfully  to  simulate  different 
engineering  problems.  This  software  has  a  specific  model  that  can 
simulate  a  range  of  different  melting  and  solidification  problems 
in  engineering,  including  casting,  melting,  crystal  growth,  and 
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Table  1 

Self-developed  numerical  model  for  different  PCM  geometries. 


Ref 

PCM 

Geometry 

Dimension 

Discretization 

methods 

Numerical  solution  methods 

Validation 

[66] 

Gallium 

Rectangular 

2D  (X,  Y) 

FDM 

Enthalpy  method-based  control  volume  approach;  the  power  law  used  for 
discretization  and  semi-implicit  method  for  pressure-linked  equations  (SIMPLE)  was 
adopted  for  the  pressure-velocity  couple  and  for  the  tridiagonal  matrix  algorithm 
(TDMA)  solver  to  solve  the  algebraic  equations 

Experimental 

[67] 

n-Octadecane 

Rectangular 

ID  (X)  2D 
(X,  Y) 

FDM 

An  enthalpy  formation-based  control  volume  approach,  the  fully  implicit  method,  and 
TDMA  were  used  to  solve  the  algebraic  equations 

[68,69] 

[70] 

Paraffin  wax 

RII-56 

PCM  bags 

In  duct 

3D  (X,  Y,  Z) 

FDM 

The  interpolating  cubic  spline  function  method  was  used  to  determine  the  specific 
heat  in  each  time  step  of  the  temperature  calculations;  the  control  volume  approach 
and  the  backward  Euler  scheme  were  chosen  for  the  temperature  discretization  in 
time 

The  power  law  used  for  discretization  and  SIMPLE  were  adopted  for  pressure-velocity 
coupling 

[71] 

[72] 

Paraffin  RT30 

Cylindrical 

2D  (R,  X) 

FVM 

Experimental 

[73] 

n-Octadecane 

Cylindrical 

2D  (R,  X) 

FDM 

Enthalpy  methods  and  the  TDMA  solver  were  used  to  solve  the  equations 

Experimental 

[74] 

4  PCM 

Cylindrical 

2D  (R,  X) 

FDM 

Fixed  control  volume  and  enthalpy  methods 

- 

[75] 

RT60 

Cylindrical 

2D  (R,  0) 

FDM 

A  control  volume-based  implicit  form  was  used.  All  equations  were  solved 
simultaneously  using  the  Gauss-Seidel  iterative  method 

Experimental 

[76] 

Paraffin  130/ 
135  Type  1 

Cylindrical 

2D  (R,  0) 

FDM 

The  enthalpy  formulation  approach,  the  control  volume  method,  and  the  alternating 
direction  scheme  were  used  to  discretize  the  basic  equations 

Experimental 

[77] 

Paraffin  wax 

Cylindrical 

2D  (R,  0) 

FDM 

The  control  volume  approach,  the  power  law  used  for  discretization,  and  the  SIMPLE 
algorithm  were  adopted  for  pressure-velocity  coupling 

[78,79] 

[80] 

n-Hexacosane 

Cylindrical 

ID  (R) 

— 

A  temperature  and  thermal  resistance  iteration  method  was  developed  to  analyze 
PCM  solidification  to  obtain  the  solutions 

Experimental 

[81] 

4PCM 

Cylindrical 

2D  (R,  X) 

FDM 

Enthalpy  methods  were  used;  the  equations  were  solved  simultaneously  using  the 
Gauss-Seidel  iterative  method 

[82] 

2PCM 

Cylindrical 

2D  (R,  X) 

FDM 

The  enthalpy  method  and  the  TDMA  were  used  to  solve  the  resulting  algebraic 
equations  solved  using  an  iteration  procedure 

[73] 

[83] 

5  PCM 

Cylindrical 

2D  (R,  X) 

FEM 

Standard  Galerkin  finite  element  method  was  used 

- 

[84] 

3  PCM 

Cylindrical 

2D  (R,  X) 

FDM 

Using  the  enthalpy  method,  the  equations  were  solved  by  adopting  the  alternating 
direction  implicit  (ADI)  method  and  the  TDMA 

[85] 

[86] 

Water 

Spherical 

ID  (R) 

FDM 

An  implicit  FDM  approach  and  a  moving-grid  scheme  were  used.  The  problem  domain 
was  divided  into  two  regions,  solid  and  liquid,  both  of  which  were  separated  by  the 
interface.  All  algebraic  equations  were  solved  simultaneously 

Experimental 

[87] 

Water,  water 

Glycol 

mixtures 

Spherical 

ID  (R) 

FDM 

The  heat  delivered  by  the  capsule  could  also  be  calculated  by  the  thermal  resistance 
circuit.  The  resulting  equation  and  associated  boundary  conditions  were  treated  using 
a  moving  grid  method 

Experimental 

FDM;  Finite  difference  method,  FEM;  Finite  elements  method,  FVM;  Finite  volume  method. 


solidification.  The  program  can  be  used  to  solve  the  phase  change 
that  occurs  at  a  single  temperature  (pure  metals)  or  over  a  range 
of  temperatures  (mixture,  alloy,  and  so  on).  The  applications  and 
limitations  of  Fluent  can  be  found  in  Ref.  [25]. 

To  begin  the  Fluent  software  simulation,  the  physical  engi¬ 
neering  problem  is  drawn  and  meshed  in  a  specific  geometric 
modeling  using  mesh  generation  tools  that  include  the  Fluent 
software  (Gambit,  workbench).  The  mesh  generation  tools  can  be 
imported  from  other  programs  like  Auto  CAD.  After  the  physical 
configuration  is  drawn  and  meshed,  the  boundary  layers  and 
zone  types  are  defined,  and  the  mesh  is  exported  to  the  Fluent 
software. 

Different  grid  sizes  and  time  steps  should  be  applied  to  the 
numerical  model  to  ensure  that  the  numerical  results  are  inde¬ 
pendent  of  the  parameters.  Small  grid  sizes  and  time  steps  are 
preferred  for  a  short  simulation  time  in  the  computer. 

2.2.  Mathematical  formulation  for  Fluent 

The  mathematical  equations  used  to  solve  the  solidification 
and  melting  models  in  Fluent  depend  on  the  enthalpy-porosity 
technique  [26-28]  and  on  the  finite  volume  methods.  In  the 
former,  the  melt  interface  is  not  tracked  explicitly.  A  quantity 
called  liquid  fraction,  which  indicates  the  fraction  of  the  cell 
volume  in  liquid  form,  is  associated  with  each  cell  in  the  domain. 
The  liquid  fraction  is  computed  at  each  iteration  based  on 
enthalpy  balance.  The  mushy  zone  is  a  region  wherein  the  liquid 
fraction  lies  between  0  and  1.  The  mushy  zone  is  modeled  as  a 
“pseudo”  porous  medium  in  which  the  porosity  decreases  from 


1  to  0  as  the  material  solidifies.  When  the  material  has  fully 
solidified  in  a  cell,  the  porosity  becomes  zero,  resulting  in  the 
drop  of  velocities  to  zero.  In  this  section,  an  overview  of  the 
solidification/melting  theory  is  given.  For  details  on  the  enthalpy- 
porosity  method,  refer  to  Voller  and  Prakash  [26]. 

The  energy  equation  is  expressed  as 

dt(pH)  +  V(pvH)  =  V(kVT)  +  S  (3) 

where  p  is  the  density  of  the  PCM,  v  is  the  fluid  velocity,  k  is  the 
thermal  conductivity,  H  is  the  enthalpy,  and  S  is  the  source  term. 
The  sensible  enthalpy  can  be  expressed  as 


h  =  /iref  +  f  cpAT, 

J  Tre  f 

(4) 

and  H  can  be  defined  as 

H  =  fi+AH, 

(5) 

where  href  is  the  reference  enthalpy  at  the  reference  temperature 
rref,  cp  is  the  specific  heat,  AH  is  the  latent  heat  content  that  may 
change  between  zero  (solid)  and  1  (liquid),  L  is  the  latent  heat  of 
the  PCM,  and  y  is  the  liquid  fraction  that  occurs  during  the  phase 
change  between  the  solid  and  liquid  state  when  the  temperature 
is  Tj  >  T  >  Ts.  Thus,  y  may  be  written  as 

y  =  A H/L,  (6) 

(  0  if  T  <  Ts 

y=\  1  if  T>Tl  (7) 

I  (T—Ts)/(T\—Ts)  if  Tl>T>Ts 
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The  source  term  Sj  in  momentum  equation,  is  defined  as 

s=c(1-y)2^r  <8> 

It  is  the  Darcy’s  law  damping  terms  (as  source  term)  that  are 
added  to  the  momentum  equation  due  to  phase  change  effect  on 
convection.  Where  e  =0.001  is  a  small  computational  constant 
used  to  avoid  division  by  zero,  the  mushy  zone  constant  C 


Pressure-Based  Segregated  Algorithm 


describes  the  kinetic  process  in  the  mushy  zone  ordinary  between 
104  and  107. 

2.3.  Fluent  solver 

The  Fluent  software  has  two  main  solvers:  the  pressure-based 
solver  and  the  density-based  coupled  solver  (DBCS).  Only  the  first 
method  can  be  used  to  simulate  the  melting  and  solidification 

Pressure-Based  Coupled  Algorithm 


Fig.  2.  Overview  of  the  pressure-based  solution  methods  [25]. 


Fig.  3.  Solution  procedure  overview  for  Fluent  [30]. 
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problems.  The  pressure-based  solver  employs  an  algorithm  that 
belongs  to  a  general  class  of  methods  called  the  projection 
method  [29]  wherein  the  constraint  of  the  mass  conservation 
(continuity)  of  the  velocity  field  is  achieved  by  solving  a  pressure 
(or  pressure  correction)  equation.  The  pressure  equation  is 
derived  from  the  continuity  and  the  momentum  equations  in 
such  a  way  that  the  velocity  field,  corrected  by  the  pressure, 
satisfies  the  continuity.  The  governing  equations  are  nonlinear 
and  coupled  to  one  another;  thus,  the  solution  process  involves 
iterations  where  the  entire  set  of  governing  equations  is  solved 
repeatedly  until  the  solution  converges.  Available  in  Fluent  are 
two  pressure-based  solver  algorithms,  namely,  a  segregated 
algorithm  and  a  coupled  algorithm  Fig.  2. 

Different  discretization  (interpolation  methods)  schemes 
are  available  for  the  convection  terms  in  Fluent,  including  the 
first-order  upwind  scheme,  power  law  scheme,  second-order 
upwind  scheme,  central-differencing  scheme,  and  the  quadratic 
upwind  interpolation  scheme.  The  first-order  upwind,  power  law, 
and  second-order  upwind  schemes  are  mostly  used  with  solidifi¬ 
cation  and  melting  problems.  The  last  two  methods  are  more 
accurate  than  the  first  methods.  In  addition,  the  interpolation 
method  of  the  face  pressure  for  the  momentum  equations  are  the 
semi-implicit  pressure-linked  equation  (SIMPLE),  the  SIMPLE- 
consistent  (SIMPLEC),  the  pressure-implicit  with  splitting  of  opera¬ 
tors  (PISO),  and  the  fractional  step  method  (FSM).  Fig.  3  shows  the 
solution  procedure  overview.  More  details  about  the  solution, 
initialization,  and  discretization  methods  can  be  found  in  Ref.  [25]. 

The  physical  properties  of  materials,  such  as  density, 
thermal  conductivity,  heat  capacity,  and  viscosity,  may  be  tem¬ 
perature-dependent  and/or  composition-dependent.  The  temperature 
dependence  is  based  on  a  polynomial,  piecewise-linear,  or 


Time,  min 


piecewise-polynomial  function.  Individual  component  properties 
are  either  defined  by  the  user  or  computed  via  kinetic  theory. 
Nevertheless,  these  physical  properties  can  be  defined  as  a 
constant  value,  a  temperature-dependent  function,  or  a  user- 
defined  function  (UDF)  that  can  be  written  in  a  specific  program¬ 
ming  language  to  define  the  temperature-dependence  of  the 
thermophysical  properties. 

Some  researchers  deduced  that  the  thermophysical  properties 
of  PCMs,  such  as  density  and  viscosity,  are  dependent  on  tem¬ 
perature  changes  and  are  determined  by  specific  correlations. 

p  =  Pi/W-Ti)  +  V  (9) 

jli  =  exp(A  +  B/T)  (10) 

Where  p\  is  the  density  of  PCM  at  the  melting  temperature  T/, 
P  is  the  thermal  expansion  coefficient,  and  A,  B  are  constant 
coefficients.  The  relationship  between  PCM  and  air  is  described  by 
a  specific  model  called  the  volume  of  fluid  (VOF).  This  model 
defines  the  PCM-air  system  with  a  moving  internal  interface,  but 
without  inter-penetration  of  the  two-phase  fluid.  Three  condi¬ 
tions  reflect  the  fluid’s  volume  fraction  in  the  computational  cells, 
denoted  as  an:  if  an  =  0,  then  the  cell  is  empty  of  the  nth  fluid;  if 
an  =  l,  then  the  cell  is  full  of  the  nth  fluid;  and  if  0<an<  1,  then 
the  cell  contains  the  interface  between  the  nth  fluid  and  one  or 
more  other  fluids. 

Shmueli  et  al.  [31]  numerically  investigated  the  melting  of 
PCM  in  a  vertical  cylinder  tube  with  a  diameter  of  3  and  4  cm, 
exposed  to  air  above,  and  insulated  at  the  bottom.  The  wall 
temperature  was  between  10  and  30  °C  above  the  melting 


Time,  min 


Fig.  4.  Melt  fraction  vs.  time  for  various  values  of  the  mushy  zone  constant  with 
different  models  [31]. 


Fig.  5.  Comparison  between  pressure  discretization  schemes  and  experimental 
results  with  different  models  [31]. 
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temperature  of  the  PCM.  They  studied  the  effect  of  various 
numerical  solution  parameters  such  as  the  pressure  velocity 
scheme  (PISO  and  SIMPLE)  and  the  pressure  discretization  regime 
(e.g.,  PRESTO  and  the  weighted  force  method).  They  analyzed  the 
influence  of  the  mushy  zone  constant  to  the  melting  process, 
as  shown  in  Fig.  4  and  Fig.  5.  They  also  reported  no  difference  in 
the  pressure  velocity  scheme  (PISO  and  SIMPLE)  as  well  as  a 
considerable  difference  in  the  result  for  pressure  discretization. 

3.  Numerical  simulation  of  PCM  in  thermal  storage 
geometries 

3.2.  Spherical  capsule  geometry 

LHTES  in  a  spherical  container  represents  an  important  case 
for  thermal  storage  to  be  used  in  different  engineering  applica¬ 
tions  such  as  in  packed-bed  storage.  Moreover,  most  commercial 
companies  produce  the  spherical  capsule  for  this  application. 

Assis  et  al.  [32]  presented  a  numerical  and  an  experimental 
investigation  of  the  melting  process  of  RT27  with  a  PCM  in  a 
spherical  geometry  filled  with  98.5%  solid  PCM  as  an  initial 
condition  of  the  simulation.  One  percent  was  left  to  account  for 
the  increase  of  the  PCM  volume  during  the  phase  change  transi¬ 
tion.  The  variable  density  was  defined  during  the  liquid  state  with 
linear  variation  in  the  “mushy”  state,  and  the  VOF  was  used  for 
the  PCM-air  system.  The  numerical  model  was  based  on  the  axial 
symmetry  of  the  physical  model  using  Fluent  6.0.  Different  design 
and  operation  parameters  were  studied,  and  the  results  indicated 
that  the  melting  process  of  the  PCM  is  affected  by  its  thermal  and 
geometrical  parameters,  including  the  Stefan  number  and  the 
shell  diameter.  Assis  et  al.  [33]  numerically  and  experimentally 
studied  the  solidification  process  of  RT27  as  a  PCM  in  a  spherical 
shell  filled  with  98.5%  liquid  PCM.  Different  shell  diameters  were 
examined,  and  a  transient  numerical  model  was  developed  using 
Fluent  6.2;  the  numerical  model  was  the  same  one  used  in  the 
melting  process  of  PCM  [32]. 

Hosseinizadeh  et  al.  [34]  numerically  studied  the  effect  of 
various  volume  fractions  of  nanopractical  copper  on  an  uncon¬ 
strained  melting  rate  in  a  spherical  container.  Different  volumes 
of  nanopractical  copper  (0,  0.02,  0.04)  per  volume  were  used,  and 
85%  of  the  PCM  was  filled.  They  used  Fluent  to  develop  an  axi- 
symmetric  numerical  model  and  added  the  Darcy  law  to  the 
momentum  equation  to  account  for  the  effect  of  phase  change  on 
convection.  The  VOF  used  for  the  PCM-air  system  and  the  power 
law  scheme  and  the  PISO  method  used  for  the  pressure-velocity 
coupling  were  adopted  to  solve  the  momentum  and  energy 
equations,  respectively.  The  PRESTO  scheme  was  used  for  the 
pressure  correction  equation.  Three  different  Stefan  numbers 
were  analyzed;  the  validation  of  the  model  was  based  on  the 
experimental  work  by  Assis  et  al.  [32].  The  results  indicated  that 
the  nanopractical  copper  caused  an  increase  in  the  thermal 
conductivity  of  the  nano-enhanced  PCMs  compared  with  the 
conventional  PCM.  Unfortunately,  the  latent  heat  of  fusions 
decreased. 

Tan  et  al.  [35]  reported  an  experimental  and  numerical  work 
on  the  constrained  melting  of  PCM  inside  a  transparent  spherical 
glass  capsule.  The  Darcy  law  was  added  to  the  momentum 
equation  to  account  for  the  effect  of  phase  change  on  convection. 
They  used  Fluent  6.2.16  for  the  numerical  simulation  and  the 
SIMPLE  method  for  solving  the  governing  equations.  The  power 
law  differencing  scheme  was  used  to  solve  the  momentum  and 
energy  equations,  whereas  the  PRESTO  scheme  was  adopted  for 
the  pressure  correction  equation.  They  reported  that  waviness 
and  excessive  melting  of  the  bottom  of  the  PCM  was  under¬ 
estimated  by  the  experimental  observation.  This  discrepancy 


was  linked  to  the  use  of  a  support  structure  to  hold  the  sphere 
that  could  have  inhibited  heat  from  reaching  the  bottom  of  the 
sphere. 

Xia  et  al.  [36]  developed  an  effective  packed-bed  thermal 
energy  storage  that  contains  a  spherical  PCM  capsule.  They 
studied  the  effect  of  the  arrangement  of  the  PCM  spheres  and 
the  encapsulation  of  the  PCM  on  the  heat  transfer  performance  of 
the  LTES  bed.  They  developed  a  2D  numerical  model  via  Fluent 
6.2  using  the  standard  k-£  model  for  turbulent  flow;  a  SIMPLE 
algorithm  was  adopted  for  pressure  and  velocity  coupling.  They 
reported  that  random  packing  was  more  favorable  than  special 
packing  for  heat  storage  and  retrieval;  both  the  material  and  the 
thickness  of  the  encapsulation  had  apparent  effects  on  the  heat 
transfer  performance  of  the  LTES  bed. 

3.2.  Square  and  rectangular  geometry 

Arasu  and  Mujumdar  [37]  presented  a  numerical  investigation 
on  the  melting  of  paraffin  wax  dispersed  with  aluminum  (A1203) 
as  a  nano-PCM  in  a  square  enclosure  heated  from  the  bottom  side 
and  from  a  vertical  side  with  different  volume  percentages 
(2%  and  5%).  They  developed  a  numerical  model  using  Fluent 
6.3.26.  The  computational  domain  was  resolved  with  a  fine  mesh 
near  the  hot  and  cold  wall  to  resolve  the  boundary  layer;  an 
increasingly  coarser  mesh  was  used  in  the  rest  of  the  domain  to 
reduce  computational  time.  They  used  a  UDF  to  define  the 
temperature-dependence  of  the  thermophysical  properties  of 
PCM,  wherein  the  first-order  upwind  difference  scheme  was  used 
to  solve  the  momentum  and  energy  equations,  and  the  PRESTO 
scheme  was  adopted  for  the  pressure  correction  equation.  They 
reported  that  the  effective  thermal  conductivity  of  a  paraffin  wax 
latent  heat  storage  medium  could  be  increased  significantly 
by  smaller  volumetric  concentrations  of  alumina  particles  in  the 
paraffin  wax. 

Pinelli  and  Piva  [38]  developed  a  2D  numerical  model  using 
Fluent  to  study  the  effects  of  natural  convection  on  the  PCM 
process  in  terms  of  temperature  distributions,  interface  displace¬ 
ment,  and  energy  storage,  and  to  enhance  the  agreement  between 
the  experimental  work  with  the  numerical  one  when  the  con¬ 
duction  dominated  is  considered  only  for  the  heat  transfer  of  the 
PCM.  The  mushy  zone  constant  for  the  model  was  approximately 
1.6*103,  and  the  cylinder  cavity  was  filled  with  n-octadecane  as 
PCM  in  the  solid  state  and  heated  from  above  with  an  electrical 
heater.  The  overall  heat  transfer  coefficient  U  considered  conduc¬ 
tion  through  the  polystyrene  layer  as  well  as  the  convection  and 
radiation  heat  exchanged  with  the  environment.  The  model  was 
validated  with  an  analytical  solution  wherein  the  conduction- 
dominated  freezing  and  numerical  methods  were  considered  in 
the  convection-conduction.  They  reported  that  the  agreement 
between  numerical  and  experimental  results  significantly  improved 
when  the  presence  of  circulation  in  the  liquid  phase,  instead  of  in 
the  conduction-dominated  process  only,  was  considered. 

Researchers  reported  that  the  heat  transfer  in  PCM  could 
be  enhanced  significantly  by  adding  nanopractical  copper  to  the 
PCM  (nanofluid).  Such  addition  improves  the  PCM’s  thermal 
conductivity.  Wu  et  al.  [39]  numerically  investigated  the  melting 
processes  of  Cu/paraffin  nanofluid  PCMs.  They  developed  a  2D 
numerical  model  using  Fluent  6.2.  The  enclosed  cavity  was  filled 
with  2.5-cm  high  PCM  and  1  cm  air;  the  top  of  the  cavity  was 
enclosed  and  taken  by  air.  The  VOF  was  adopted  to  describe  the 
relationship  between  the  PCM  and  air  with  a  moving  interface, 
but  without  inter-penetration  of  the  two  media.  The  QUICK 
differencing  scheme  was  used  to  solve  the  momentum  and  energy 
equations,  whereas  the  PRESTO  scheme  was  adopted  for  the 
pressure  correction  equation.  They  reported  a  13.1%  decrease  in 
the  paraffin’s  melting  time. 
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Li  et  al.  [40]  numerically  studied  the  melting  process  of  a  PCM 
in  concentric  horizontal  annuli  of  a  square  external  shell.  A  2D 
numerical  model  was  adopted  using  Fluent  6.0.  They  reported 
that  the  top  part  of  the  PCM  melts  faster  than  the  bottom  part  as  a 
result  of  the  convection-dominated  mode,  which  accelerates  the 
melting  of  the  top  part  of  the  annulus. 

Khodadadi  and  Hosseinizadeh  [41]  reported  that  the  disper¬ 
sion  of  nanoparticles  in  PCMs  enhanced  the  material’s  thermal 
conductivity  in  contrast  to  the  base  material.  They  used  a  2D 
numerical  model  within  the  commercial  code  of  Fluent  6.2.1.  The 
SIMPLE  algorithm  and  the  QUICK  differencing  scheme  were  used 
to  solve  the  momentum  and  energy  equations,  whereas  the 
PRESTO  scheme  was  adopted  for  the  pressure  correction  equation. 
The  Darcy  Law  damping  terms  were  added  to  the  momentum 
equation.  The  results  indicated  an  increase  in  the  thermal  con¬ 
ductivity  of  the  nano-enhanced  PCMs  in  contrast  to  the  conven¬ 
tional  PCM.  Unfortunately,  the  latent  heat  of  fusion  decreased. 

3.3.  Cylindrical  capsule  geometry 

Cylindrical  geometries  are  considered  most  promising  for 
devices  for  commercial  heat  exchangers,  such  as  the  double  pipe 
heat  exchanger  and  shell  and  the  tube  heat  exchanger,  because  of 
their  high-efficiency  in  a  minimum  volume.  Most  researchers  use 
this  geometry  for  LHTES,  filling  the  PCM  in  the  tube  side  or  in  the 
annulus  (shell).  Commercial  products  of  PCMs  are  delivered  in 
these  geometries. 

Darzi  et  al.  [42]  used  Fluent  6.2  to  simulate  the  melting 
process  of  the  n-octadecane  as  a  PCM  in  the  concentric  and 
eccentric  double  pipe  heat  exchanger.  They  studied  the  effect  of 
the  internal  tube  position  on  the  melting  rate.  They  developed  a 
2D  numerical  model  (R,  0)  and  adapted  low  power  to  solve  the 
momentum  and  energy  equations.  The  SIMPLE  method  was  used 
for  pressure-velocity  coupling,  and  the  PRESTO  scheme  was  used 
for  the  pressure  correction  equation.  The  results  indicated  that 
the  melting  rate  was  the  same  in  the  early  stage  of  melting,  and 
then  decreased  in  the  concentric  model. 

Gou  and  Zhang  [43]  reported  that  the  solidification  process 
of  KN03-NaN03  as  PCM  in  solar  power  generation  using  the 
direct  steam  technology  is  controlled  by  pure  conduction.  They 
neglected  the  influences  of  natural  convection  when  comparing  it 
with  heat  conduction.  They  developed  a  2D  numerical  model  via 
Fluent  6.2  and  studied  the  storage  with  and  without  foil  based  on 
the  same  tube  diameters  and  PCM.  The  results  indicated  how  the 
heat  transfer  and  discharge  time  were  affected  by  the  changes  in 
the  geometry  of  the  aluminum  foil  and  in  the  tube  radius.  They 
reported  that  the  discharge  process  significantly  accelerated  with 
the  addition  of  the  aluminum  foil. 

Colella  et  al.  [2]  designed  a  medium-scale  LF1TES  unit  for 
district  heating  systems.  The  shell-and-tube  LF1TES  that  included 
RT100  as  a  PCM  in  the  shell  side  was  used  to  transfer  heat  from 
the  district  to  the  heating  networks  of  the  building;  the  heat 
transfer  fluid  was  water.  The  thermal  storage  was  charged  during 
the  night  using  CLIP  plants  and  then  discharged  during  the  day  at 
the  peak  of  the  thermal  request.  A  2D  numerical  model  (R,  X)  was 
developed  using  the  Fluent  simulation  software  to  study  the 
thermal  behavior  of  the  PCM  and  the  HTF.  Grid  meshing  of  the 
computational  model  was  refined  systematically.  The  SIMPLE 
algorithm  was  adopted  to  solve  the  pressure-velocity  coupling. 
The  second-order  upwind  scheme  was  used  for  the  convective 
fluxes.  A  naturally  expanded  graphite  matrix  with  15%  volume 
was  added  to  the  paraffin  to  enhance  the  thermal  conductivity  of 
the  composite  (paraffin-graphite)  to  achieve  the  heat  fluxes 
required  for  the  medium  scale  application.  Different  operation 
strategies  and  parameters  were  studied,  including  the  time-wise 
variation  of  the  liquid  fraction,  PCM  temperature,  and  the  HTF 


outlet  temperature.  Based  on  their  analysis,  the  best  options  for 
energy  storage  density  were  characterized  by  latent  storage  units 
installed  on  the  heating  network  of  the  building  rather  than  on 
the  primary  district  heating  network. 

Buyruk  et  al.  [44]  studied  the  solidification  process  of  water 
around  different  staggered  and  inline  cylinders  placed  in  a 
rectangular  ice  storage  tank.  They  developed  a  2D  numerical 
model  using  Fluent  to  depict  the  temperature  distributions  and 
ice  formation  in  the  tank.  They  adopted  Darcy’s  law  and  the 
Kozeny-Carman  equation  to  model  the  flow  and  permeability  of 
the  mushy  zone,  respectively.  Based  on  a  fixed  ice  storage  tank 
volume,  the  solidified  area  of  four  inline  cylinders  was  larger  than 
that  for  six  inline  cylinders. 


4.  Numerical  simulation  of  the  PCM  applications 

4.1.  Numerical  simulation  of  PCM  for  electronic  devices  applications 

Shatikan  et  al.  [45]  numerically  presented  the  melting  process 
of  a  PCM  with  internal  fins  exposed  to  air  from  the  top  and  heated 
from  the  horizontal  base.  They  investigated  different  fin  dimen¬ 
sions  with  a  constant  ratio  between  the  fin  and  the  PCM  thick¬ 
ness.  They  used  Fluent  6.0  to  develop  2D  and  3D  numerical 
models  and  adopted  the  VOF  model  to  describe  the  relationship 
between  the  air  and  the  PCM;  the  SIMPLE  algorithm  was  used  for 
the  pressure-velocity  coupling.  The  2D  model  was  considered 
because  the  3D  model  was  time  consuming.  They  reported  that 
the  results  for  the  2D  and  3D  models  were  equal  in  some 
dimensions,  and  the  transient  phase  change  process  of  the  PCM 
depended  on  the  operation  and  design  parameters  of  the  system. 
Such  parameters  were  related  to  the  temperature  difference 
between  the  base  and  the  mean  melting  temperature,  as  well  as 
to  the  thickness  and  height  of  the  fins.  Shatikian  et  al.  [10] 
analyzed  the  same  model  with  constant  heat  flux  applied  from 
the  horizontal  base. 

Wang  and  Yang  [11]  numerically  studied  the  cooling  technol¬ 
ogy  of  portable  handheld  electronic  devices  using  PCM.  They 
developed  half  computational  grids  with  the  3D  numerical  model 
using  Fluent  12  to  investigate  the  effect  of  a  different  amount  of 
fins  and  various  power’s  heating  level.  They  adopted  the  VOF 
model  to  describe  the  relationship  between  the  air  and  the  PCM. 
In  the  VOF  model,  the  Geo-Reconstruct,  which  is  the  volume 
fraction  discretization  scheme,  was  used  to  calculate  the  transient 
moving  boundary.  If  the  transient  temperature  of  the  calculation 
cell  in  the  PCM  domain  was  equal  to  or  higher  than  the  melting 
temperature,  the  continuity  and  momentum  equations  were 
calculated  in  the  system;  the  SIMPLE  algorithm  was  used  for  the 
pressure-velocity  coupling.  They  studied  the  orientation  of  PCM 
with  the  heat  sink;  the  model  was  validated  using  the  method  by 
Fok  et  al.  [46].  The  numerical  results  indicated  that  the  orienta¬ 
tion  had  insufficient  influence  on  the  transient  thermal  perfor¬ 
mance  of  the  hybrid  cooling  system,  and  that  the  transient  surface 
temperature  caused  a  discrepancy  of  around  8.3%-10.7%.  They 
reported  that  PCM  in  aluminum  heat  sink  with  six  fins  could 
provide  stable  temperature  control  and  good  cooling  ability. 

Yang  and  Wang  [12]  numerically  investigated  the  cooling 
technique  using  a  hybrid  PCM-based  heat  sink  within  a  closed 
system  to  absorb  constant  and  uniform  volumetric  heat  generated 
from  portable  handheld  electronic  devices.  They  tested  different 
computation  grids  and  time  steps  to  validate  the  model  with 
the  experimental  work;  various  power  levels  (2-4  W),  different 
orientations  (vertical,  horizontal,  slanted),  as  well  as  charging 
and  discharging  times  were  conducted  in  the  simulation  model. 
The  melting  rate  increased  when  a  higher  power  level  was 
transferred  to  the  system.  They  reported  that  the  maximum 
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temperature  of  the  hybrid  system  using  PCM  could  be  well 
controlled  under  320  K  if  the  melting  ratio  was  under  0.9  or  if 
the  melting  time  was  shorter  than  6801.2  s. 

Ye  et  al.  [47]  figured  out  the  correlation  between  influential 
variables,  such  as  the  liquid  fraction  and  thermal  storage  time,  the 
transient  heat  flux  and  melting  period,  and  the  solid  fraction  and 
the  solidification  time  by  numerically  studying  the  heat  transfer 
and  fluid  flow  in  a  plate  fin  filled  with  PCM  for  rapid  heat 
absorbed/released  techniques.  They  developed  a  2D  numerical 
model  through  Fluent.  The  cavity  was  filled  with  85%  PCM  and 
15%  air  to  consider  thermal  expansion  during  the  solid-liquid 
phase  change;  finely  structured  meshing  grids  were  installed  near 
the  heating  and  cooling  plate  wall  to  resolve  the  boundary  layers. 
The  VOF  model  was  adopted  to  describe  the  relationship  between 
the  air  and  the  PCM.  The  UDF  was  written  in  C  language  to 
account  for  the  temperature  dependence  of  the  thermophysical 
properties.  The  SIMPLE  algorithm  was  considered  for  the  pressure- 
velocity  coupling,  and  different  time  steps  were  used  for  the 
energy  storage  and  energy  release.  Exploring  via  a  3D  numerical 
model,  they  reported  no  difference  in  the  results  between  the  2D 
and  3D  models.  The  correlations  they  obtained  will  be  useful  for 
future  component  design  and  system  optimization. 

Ye  et  al.  [48]  performed  further  studies  on  the  effect  of  PCM 
cavity  volume  fraction  on  the  heat  transfer  behavior  and  fluid 
flow  in  LF1TES.  Different  parameters  of  LHTES  performance,  such 
as  the  volume  expansion  ratio,  the  time  of  complete  thermal 
storage,  heat  flux,  liquid  fraction,  and  velocity  and  temperature 
fields,  were  investigated.  They  used  the  same  model  in  Ref.  [47], 
except  that  the  time  step  was  chosen  in  relation  to  the  grid  size 
via  the  Courant-Friedriches-Lewy  number.  They  reported  that 
the  volume  expansion  ratio  decreased  as  cavity  volume  fraction 
x  increased.  By  contrast,  the  time  for  complete  energy  storage 
increased. 

Kandasamy  et  al.  [49]  numerically  and  experimentally  inves¬ 
tigated  the  use  of  a  PCM-based  heat  sink  for  quad  flat  package 
(QAP)  electronic  devices.  They  developed  one  quarter  of  the 
physical  model  as  well  as  a  3D  numerical  model  using  the 
commercial  Fluent  6.2.  The  VOF  was  used  to  describe  the  relation¬ 
ship  between  the  air  and  the  PCM.  The  PISO  scheme  was  used  for 
the  pressure-velocity  coupling.  The  cooling  performance  of  the 
PCM-based  heat  sinks  improved,  compared  with  the  case  without 
PCM.  The  full  melting  of  the  PCM  and  the  melting  rate  also 
increased  as  the  input  power  increased. 

Sabbah  et  al.  [50]  presented  the  thermal  and  hydraulic 
behavior  of  a  micro-channel  heat  exchanger  with  MCPCM  water 
slurry  for  heat  dissipation  of  a  high-power  electronic  device.  They 
developed  a  3D  numerical  model  using  Fluent  6.2,  and  the  SIMPLE 
algorithm  was  used  for  velocity-pressure  coupling.  The  heat 
transfer  in  both  the  cooling  fluid  and  the  metal  heat  sink  on  the 
overall  performance  of  the  system  was  considered.  The  results 
showed  a  significant  increase  in  the  heat  transfer  coefficient  that 
depended  on  the  channel  inlet  and  outlet  and  on  the  selected 
MCPCM  melting  temperatures.  Lower  and  more  uniform  tem¬ 
peratures  throughout  the  electronic  device  could  be  achieved 
using  less  pumping  power  compared  with  the  case  using  only 
water  as  the  cooling  fluid. 

Flosseinizadeh  et  al.  [51]  numerically  and  experimentally 
presented  the  thermal  management  of  electronic  products  using 
PCM.  They  compared  the  effect  of  the  PCM  cooling  technology  in 
contrast  to  the  normal  ones.  Different  design  parameters  were 
studied  to  determine  their  effects  on  the  temperature  of  the  heat 
sink  and  on  the  melting  of  PCM  inside  the  heat  sink.  85%  of  the 
heat  sink  height  was  filled  with  PCM,  and  15%  of  the  height  was 
modeled  in  the  air  region  for  the  solid-to-liquid  conversion  in  the 
PCM  expansion.  A  2D  numerical  model  was  adopted  via  Fluent. 
Density  was  considered  constant  in  the  solid  state,  whereas  liquid 


density  was  a  function  of  the  temperature  and  the  thermal 
expansion  coefficient.  The  VOF  was  adopted  to  describe  the 
PCM-air  system.  The  PISO  algorithm  was  used  for  the  pressure- 
velocity  coupling.  They  set  the  time  step  at  0.0001  s  at  the  early 
stage  of  simulation  and  changed  it  to  0.001  s  after  3  min  of 
simulation.  The  increased  fin  thickness  showed  only  a  slight 
improvement  in  thermal  performance,  whereas  an  increase  in 
the  number  of  fins  and  fin  height  resulted  in  an  appreciable 
increase  in  the  overall  thermal  performance. 


4.2.  Numerical  simulation  of  PCM  for  HVAC  equipment 

Tay  et  al.  [52]  modeled  different  tube  configurations  in  a  shell 
and  tube  heat  exchanger  with  PCM  in  the  shell  side.  They 
developed  a  3D  CFD  model  using  ANSYS  code  to  analyze  the 
transient  heat  transfer  during  the  phase  change  process.  Three 
grid  resolutions  were  evaluated  for  the  tube,  HTF,  and  PCM. 
ANSYS  CFX  12.1  was  used  to  complicate  the  physical  properties 
of  the  mixture.  Three  domains  were  created  using  CFX.PRE 
version  12.1,  namely,  liquid  for  HTF,  liquid  for  PCM,  and  solid 
for  the  tube.  An  experimental  apparatus  was  developed  to 
compare  the  predicted  simulation  results  of  both  freezing  and 
melting  processes.  They  found  that  the  results  of  the  CFD  melting 
model  were  generally  longer  than  the  experimental  results.  The 
thermal  behavior  of  the  freezing  and  melting  processes  of  PCM 
was  different  from  those  in  the  experimental  results  when  the 
natural  convection  was  ignored. 

Antony  and  Velraj  [53]  studied  the  heat  transfer  and  fluid  flow 
through  a  tube  in  a  regenerative  PCM  heat  exchanger  (shell  and 
tube)  module  incorporated  into  a  ventilation  system.  The  PCM 
was  encapsulated  in  the  shell  side  of  the  module  for  a  free  cooling 
system.  Fluent  software  was  used  to  analyze  the  heat  transfer  and 
flow  through  the  model,  which  consisted  of  an  air  flow  passenger 
and  two  air  spacers  on  top  and  at  the  bottom.  PISO  algorithm  was 
adopted  for  the  pressure-velocity  coupling.  The  second-order 
upwind  scheme  was  used  for  the  approximation  of  the  convective 
terms.  The  k-£  model  was  used  for  the  turbulence  modeling  of 
the  heat  transfer  fluid.  This  module  was  validated  experimentally. 
The  results  from  the  experiment  and  the  calculated  results  were 
compared. 

Xia  and  Zhang  [54]  prepared  an  acetamide/expanded  graphite 
composite  (AC/EG)  with  10%  (mass  fraction)  as  PCM  to  be  used  in 
active  solar  applications,  such  as  solar-driven  solid/liquid  desic¬ 
cant  dehumidification  systems  and  solar-driven  adsorption/ 
absorption  refrigeration  systems  that  require  a  higher  heat  source 
temperature  of  over  60  °C.  They  measured  the  thermal  properties, 
including  melting/freezing  points,  thermal  conductivities,  and  the 
heat  of  fusion,  using  the  transient  hot-wire  method  and  a 
differential  scanning  calorimeter.  Two-thirds  of  the  test  model 
was  filled  with  solid  PCM  for  the  thermal  expansion  during  the 
solid-liquid  phase  change.  They  developed  a  2D  numerical  model 
using  Fluent  6.2  to  study  the  release  and  absorption  of  heat 
in  LHTES.  The  numerical  model  was  validated  experimentally. 
The  thermal  conductivity  of  the  composite  AC/EG  (2.61)  was  six 
times  that  of  pure  metal;  0.043  for  solid  at  60  °C  and  0.2  for  liquid 
at  90  °C.  The  phase  change  temperature  shifted  and  the  latent 
heat  decreased  compared  with  that  of  pure  AC.  The  heat  storage 
and  retrieval  durations  of  the  LTES  unit  with  the  AC/EG  composite 
showed  45%  and  78%  reductions,  respectively,  compared  with 
those  for  pure  AC. 

New  internal  fins  inside  spherical  and  cylindrical  PCM  encap¬ 
sulation  in  the  HVAC  system  for  buildings  were  numerically 
investigated  by  Siva  et  al.  [55].  They  developed  a  numerical 
model  using  Fluent  6.3  to  study  the  charging  and  discharging 
process  of  PCMs  for  the  two  configurations.  Incorporating  fins  for 
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different  configurations  reduced  the  total  solidification  time  by 
65%-72%. 

4.3.  Numerical  simulation  of  PCM  for  building  applications 

LHTES  is  used  extensively  in  building  structures  to  remove 
thermal  load  and  to  improve  the  thermal  comfort  of  occupants. 
PCM  can  be  found  in  different  items  of  the  building  structure, 
such  as  the  walls,  roof,  floor,  ceiling,  window  shutters,  plaster¬ 
board,  and  tiles.  Most  researchers  mentioned  in  the  present  study 
used  CFD  to  simulate  the  process,  reporting  a  good  agreement 
between  the  simulation  and  the  experimental  results. 

Silva  et  al.  [56]  incorporated  a  PCM  into  a  typical  clay  brick 
masonry  enclosure  wall  in  Portugal  as  a  passive  technology  to 
enhance  the  building’s  energy  efficiency.  They  experimentally  and 
numerically  evaluated  how  the  PCM  reduced  internal  tempera¬ 
ture  fluctuations  and  increased  the  time  delay  between  internal 
and  external  conditions.  They  generated  the  mesh  grid  in  CFD  and 
developed  a  2D  numerical  model  in  ANSYS  Fluent  12.  The  SIMPLE 
algorithm  was  used  to  solve  the  governing  equations.  The  second- 
order  upwind  scheme  was  used  to  solve  the  pressure,  momen¬ 
tum,  and  energy  equations.  They  reported  that  incorporating  the 
PCM  into  the  wall  contributed  to  the  attenuation  of  the  indoor 
space  temperature  swing,  which  reduced  from  10  to  5  °C,  and  to 
the  thermal  amplitude.  The  time  delay  also  increased  by  approxi¬ 
mately  3  h. 

Joulin  et  al.  [6]  analyzed  the  energy  conservation  of  a  building 
with  PCM  conditioned  in  a  parallel  epipedic  polyfin  envelope  to 
be  used  in  passive  solar  walls.  They  developed  a  ID  Fortran 
numerical  model  and  a  2D  model  using  commercial  Fluent  to 
study  the  thermal  behavior  of  the  PCM  to  validate  the  experi¬ 
mental  work.  The  self-developed  numerical  model  was  based  on 
the  enthalpy  formulation  and  on  the  fully  implicit  finite  differ¬ 
ence  solution  method  used  for  discretization.  The  equations  were 
solved  using  TDMA.  The  liquid  fraction  was  updated  in  each 
iteration  until  convergence  was  achieved.  For  the  commercial 
Fluent,  the  SIMPLE  method  was  used  to  solve  the  governing 
equations.  The  power  law  differencing  scheme  was  used  to  solve 
the  momentum  and  energy  equations,  whereas  the  PRESTO 
scheme  was  adopted  for  the  pressure  correction  equation.  They 
reported  that  the  ID  and  2D  numerical  simulations  did  not  satisfy 
because  the  results  were  very  far  from  the  experimental  ones:  the 
codes  did  not  represent  the  supercooling  phenomenon.  Numerical 
modeling  showed  that  the  supercooling  phenomenon  must  be 
taken  into  account  correctly  to  predict  the  PCM  thermal  behavior. 

Susman  et  al.  [57]  constructed  four  150  mm2  prototype  PCM 
sails  from  paraffin/low-density-polyethylene  composites  installed 
below  the  ceiling  of  an  occupied  office  space  in  London  in  the 
summer.  During  daytime,  the  modules  were  hung  on  an  internal 
rig  2.5  m  above  the  ground;  after  8  p.m.,  the  modules  were 
transferred  to  the  outside  rig  to  discharge  the  energy,  and  were 
moved  to  the  inside  rig  at  7  a.m.  They  developed  a  semi-empirical 
model  of  the  modules  in  Fluent  using  an  enthalpy-porosity 
formulation  to  model  phase  change.  The  modules  were  validated 
well  using  the  temperature  measurements,  including  a  notable 
divergence  when  the  maximum  liquid  fraction  was  reached. 

4.4.  Numerical  simulation  for  other  applications 

Bo  et  al.  [58]  introduced  a  metal  PCM  energy  storage  in  a 
high-temperature  heat  pipe  steam  boiler  used  for  direct  steam 
generation  in  a  solar  thermal  power  generation  system.  The  Al-Si 
eutectic  alloy  contained  12.07%  Si  used  as  a  PCM  at  a  melting 
temperature  of  577  °C  The  temperature  distributions  of  the  heat 
reservoir  tank  and  the  heat  storage  medium  under  high  solar  heat 
flux  were  analyzed  using  a  2D  numerical  model  via  Fluent.  The 


newly  developed  heat  storage  boiler  could  withstand  a  focused 
solar  energy  flux  of  400  kW/m2  and  was  compatible  with  the  heat 
storage  medium. 

Tan  et  al.  [16]  used  water  as  a  PCM  to  recover  and  store  cold 
energy  from  a  cryogenic  gas  in  a  cryogenic  cold  energy  system. 
They  developed  a  2D  numerical  model  using  Fluent  to  analyze  the 
effect  of  the  thermal  boundary  and  thermal  resistance  on  the 
freezing  process  in  an  LHTES  system.  Their  numerical  model  was 
based  on  half  the  domain  because  the  geometrical  model  was 
considered  symmetric.  The  standard  /<-{  model  was  adopted  to 
calculate  the  internal  forced  convection  of  the  gaseous  coolant. 
The  CFD  model  was  validated  experimentally  and  showed  good 
agreement  with  the  experimental  results.  They  reported  that  the 
dimensionless  numbers,  such  as  the  Biot  number  and  the  Stefan 
number  of  the  PCM,  and  the  Stanton  number  of  the  coolant 
flowing  in  the  tube  had  a  noticeable  influence  on  the  PCM 
freezing  characteristics.  A  larger  Biot  number  was  beneficial  to 
the  heat  transfer  between  the  PCM  and  the  coolant,  promoting  a 
higher  freezing  rate.  A  higher  Stefan  number  appeared  to  result  in 
larger  freezing  rates  in  a  fixed  axial  position.  Moreover,  the  frozen 
layer  slope  along  the  tube  was  steeper  at  larger  Stanton  numbers. 

Xiaohong  et  al.  [59]  numerically  studied  a  PCM  at  high  melting 
temperatures  (11 20  K)  in  a  heat  pipe  receiver  in  an  advanced 
solar  dynamic  power  generation  system.  A  2D  numerical  model 
was  adopted  via  Fluent  6.2  to  simulate  the  heat  transfer  of  the 
PCM  canister  of  the  heat  pipe  receiver.  UDF  was  used  in  the 
simulation.  The  numerical  results  were  compared  with  NASA’s 
numerical  results.  The  temperature  of  the  heat  pipe  wall  and 
the  outer  wall  of  the  PCM  canister  in  the  numerical  model  and  in 
the  NASA  simulation  differed  only  slightly.  These  results  can  be 
used  to  guide  the  design  and  optimization  of  PCM  canisters  for 
heat  pipe  receivers. 

Koizumi  and  Jin  [60]  proposed  a  new  compact  slab  container  with 
an  arc  outer  configuration  to  promote  direct  contact  melting.  They 
studied  two  configurations  of  the  close  contact  melting  process, 
namely,  flat  slab  container  and  curved  slab  container.  The  2D 
numerical  model  was  performed  via  Fluent  6.3.  The  VOF  model  was 
adopted  to  describe  the  PCM-air  system.  A  first-order  accurate 
upwind  difference  scheme  was  adopted  for  the  convection  terms 
because  the  convection  velocity  in  the  liquid  PCM  was  extremely 
small  (within  several  mm/s).  The  simulated  results  quantitatively 
elucidated  the  experimental  melting  process  from  beginning  to  end. 

MacPhee  and  Dincer  [61]  numerically  investigated  and  ther¬ 
modynamically  analyzed  the  melting  and  solidification  process  of 
a  de-ionized  water  as  PCM  in  a  spherical  capsule  geometry.  They 
developed  a  3D  model  using  Fluent  6.0  and  studied  the  effect  of 
the  HTF  inlet  temperature  and  the  flow  rate.  The  numerical  model 
was  validated  and  was  found  to  be  in  good  agreement  with  the 
experimental  data.  The  variations  of  the  incoming  HTF  tempera¬ 
ture  had  a  much  larger  effect  on  the  charging  times  and  on  the 
energy  and  exergy  efficiencies  compared  with  the  changing  of  the 
HTF  flow  rate. 

Reddy  [62]  presented  a  solar-integrated  collector  storage 
water  heating  system  with  PCM  to  optimize  the  solar  gain  and 
heat  loss  characteristics  through  a  2D  numerical  model  using 
commercial  Fluent  and  UDF  to  express  the  variable  heat  flux 
condition  on  the  absorber  plate.  Different  numbers  of  fins  were 
attached  to  the  PCM  water  solar  system  to  enhance  the  heat 
transfer.  The  latent  heat  storage  with  nine  fins  was  optimal  at 
maximum  water  temperature  and  had  minimum  heat  loss  to  the 
surroundings. 

Rao  et  al.  [63]  discussed  the  use  of  PCM  with  an  aging 
commercial  LiFeP04  power  battery  to  study  the  thermal  behavior 
of  the  PCM  and  the  cell  for  thermal  energy  management. 
They  developed  a  3D  computation  model  for  a  single  cell  and  a 
battery  package  using  the  commercial  computational  Fluent.  They 
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reported  that  the  design  for  battery  thermal  management  should 
include  a  proper  selection  of  the  thermal  conductivities  for  the 
PCM  and  the  cell  because  the  thermal  conductivity  increased, 
whereas  the  latent  heat  of  fusion  decreased. 

Darkwa  and  Su  [64]  studied  the  effect  of  particle  distribution  on 
the  thermal  performance  of  micro-encapsulated  PCMs  in  different 
configurations  of  a  composite  high  conductivity  laminated  MCPCM 
board.  They  developed  a  3D  numerical  model  using  Fluent  6.3.  The 
SIMPLE  algorithm  was  used  for  the  pressure-velocity  coupling,  and 
the  second-order  upwind  scheme  was  adopted  for  both  momentum 
and  energy  equations.  The  PRESTO  pressure  discretization  scheme 
was  used  for  the  pressure  correction  equation.  The  thermal  response 
times  for  the  rectangular  and  triangular  geometries  were  approxi¬ 
mately  half  of  that  for  the  pyramidal  geometry  during  the  cooling 
and  heating  processes  of  the  board. 

MacPhee  et  al.  [65]  numerically  studied  the  effect  of  different 
capsule  geometries,  HTF  inlet  temperatures,  and  HTF  flow  rates 
on  energy  and  exergy  efficiencies  for  the  solidification  process  in 
encapsulated  ice  thermal  energy  storage  (EITES).  They  used  Fluent 
6.3  to  simulate  different  geometries  (slab,  cylindrical,  spherical), 
three  HFT  flow  rates,  and  three  HTF  inlet  temperatures  to  achieve 
the  highest  efficient  charge  method  for  the  thermal  storage.  They 
recommended  the  exergy  efficiency  methods  to  study  system 
performance,  stating  that  these  methods  provide  better  insight 
into  system  losses.  They  also  advised  EITES  designers  to  increase 
both  the  flow  rate  and  inlet  HTF  temperature  to  achieve  full 
system  charging  in  an  acceptable  amount  of  time. 

5.  Conclusion 

The  sustainable  thermal  energy  storages  for  different  engineering 
applications,  such  as  building  structures  and  HVAC  equipment,  are 
indispensable  to  reducing  greenhouse-gas  emission.  LHTES  is  an 
important  component  of  thermal  solar  engineering  systems,  playing 
a  key  role  in  improving  the  energy  efficiency  of  these  systems. 
Thermal  energy  storage  systems,  especially  LHTES,  have  gained 
widespread  attention  in  relation  to  global  environmental  problems 
and  energy-efficiency  improvement.  The  following  conclusions  can  be 
summarized  for  CFD  applications  in  latent  heat  thermal  storage. 

1.  The  numerical  solution  for  the  PCM  phenomena  is  more 
accurate  than  the  analytical  solution  and  can  be  used  in 
different  engineering  conditions. 

2.  The  CFD  software  is  used  successively  to  simulate  the 
application  of  PCM  in  different  engineering  applications, 
such  as  electronic  cooling  technology,  building  thermal 
storages,  and  HVAC. 

3.  Different  CFD  software  are  used  in  PCM  engineering,  including 
ANSYS  Fluent,  Comsol  Multiphysics,  and  Star-CCM+.  The 
most  commonly  used  software  is  ANSYS  Fluent. 

4.  The  numerical  results  of  the  2D  numerical  model  are  generally 
the  same  as  those  of  the  3D  numerical  model.  Researchers 
reported  that  the  time  required  for  simulation  was  reduced. 

5.  The  application  of  CFD  in  designing  PCM  thermal  storages 
is  a  feasible  method  because  of  the  highly  accurate  results 
of  the  experimental  work.  CFD  also  delivers  optimization 
tools  to  help  users  achieve  maximum  efficiency  while 
saving  time  and  money. 
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